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ABSTRACT 

Cosmic rays (CRs) play a decisive role within our own Galaxy. They provide partial pressure support against gravity, they trace 
past energetic events such as supernovae, and they reveal the underlying structure of the baryonic matter distribution through their 
interactions. To study the impact of CRs on galaxy and cosmic structure formation and evolution, we develop an approximative 
framework for treating dynamical and radiative effects of CRs in cosmological simulations. Our guiding principle is to try to find a 
balance between capturing as many physical properties of CR populations as possible while at the same time requiring as little extra 
computational resources as possible.We approximate the CR spectrum of each fluid element by a single power-law, with spatially 
and temporally varying normalisation, low-energy cut-off, and spectral index. Principles of conservation of particle number, energy, 
and pressure are then used to derive evolution equations for the basic variables describing the CR spectrum, both due to adiabatic 
and non-adiabatic processes. The processes considered include compression and rarefaction, CR injection via shocks in supernova 
remnants, injection in structure formation shock waves, in-situ re-acceleration of CRs, CR spatial diffusion, CR energy losses due to 
Coulomb interactions, ionisation losses, Bremsstrahlung losses, and, finally, hadronic interactions with the background gas, including 
the associated y-ray and radio emission due to subsequent pion decay. We show that the formalism reproduces CR energy densities, 
pressure, and cooling rates with an accuracy of ~ 10% in steady state conditions where CR injection balances cooling. It is therefore 
a promising formulation to allow simulations where CR physics is included. Finally, we briefly discuss how the formalism can be 
included in Lagrangian simulation methods such as the smoothed particle hydrodynamics technique. Our framework is therefore well 
suited to be included into numerical simulation schemes of galaxy and structure formation. 

Key words. Cosmic rays - galaxies: ISM - intergalactic medium - galaxies: cluster: general - acceleration of particles - radiation 
mechanisms: non-thermal 



-""j 1. Introduction populations are an important galactic storage for energy from 
. , supernova explosions, and thereby help to maintain dynamical 
1.1. Motivation feedback of star formation for periods longer than thermal gas 
The interstellar medium (ISM) of galaxies is highly com- physics alone would permit. The spectral distribution of CRs is 
plex, with an energy budget composed both of thermal and "^'^'^h broader than that of thermal populations, which has to be 
non-thermal components. Each of the non-thermal components ^^^en into account in estimating their dynamical properties. The 
which are magnetic fields and cosmic rays (CRs) are known to dynamical important part of CR spectral distributions in mo- 
contribute roughly as much energy and pressure to the ISM as "^^"t'^"^ ^P^'^^^ ^^^^i^y encompasses a few orders of magnitude, 
the thermal gas does, at least in our own Galaxy. whereas thermal distributions appear nearly mono-energetic on 

Magnetic fields permeate the ISM and seem to have ordered logarithmic scales (see Fig. 1). 
and turbulent field components. Very likely they play an impor- Numerical simulations and semi-analytical models of galaxy 
tant role in moderating molecular cloud collapse and star forma- and large-scale structure formation neglected the effects of CRs 
tion. They are the ISM ingredient which couples otherwise dy- and magnetic fields so far, despite their dynamical importance, 
namically independent ingredients like the ISM plasma, the CR Although there have been attempts to equip SPH galaxy forma- 
gas, and (charged) dust particles into a single, however complex tion codes with magnetic field descriptions on the MHD level 
fluid. (Dolag et al. 1999), a fully dynamical treatment of the CR com- 
CRs behave quite differently compared to the thermal gas. ponent has not yet been attempted due to the very complex CR 
Their equation of state is softer, they are able to travel actively physics involved. There have been pioneering efforts to furnish 
over macroscopic distances, and their energy loss time-scales grid-based MHD codes with a diffusive CR component in order 
are typically larger than the thermal ones. Furthermore, roughly to study isolated effects like Parker instabilities (Kuwabara et al. 
half of their energy losses are due to Coulomb and ionisation 2004; Hanasz & Lesch 2003) or the injection of CRs into the 
interactions and thereby heat the thermal gas. Therefore, CR wider IGM (Miniati et al. 2001; Ryu et al. 2003; Ryu & Kang 
2003, 2004) and ISM (Snodin et al. 2006). However, these codes 
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cal setting due to the missing adaptive resolution capability, and 
the lack of a treatment of cosmological components such as dark 
matter and stellar populations. There have also been numerical 
implementations of discretised CR energy spectra on top of grid- 
based cosmological simulations (Miniati 2001), but these imple- 
mentations neglected the hydrodynamic pressure of the CR com- 
ponent. In addition, the amount of computer resources required 
for these models in terms of memory and CPU-time is substan- 
tial. This renders this approach unattractive for self-consistent 
simulations of galaxy formation and inhibits the inclusion of 
CRs into cosmological simulations of structure formation, where 
it is not clear ab initio if CRs are crucial or not. 

An accurate description of CRs should follow the evolution 
of the spectral energy distribution of CRs as a function of time 
and space, and keep track of their dynamical, non-linear cou- 
pling with the hydrodynamics. In order to allow the inclusion 
of CRs and their effects into numerical simulations and semi- 
analytic descriptions of galaxy formation, we develop a simpli- 
fied description of the CR dynamics and physics. We seek a com- 
promise between two opposite requirements, namely: (i) to cap- 
ture as many physical properties and peculiarities of CR popula- 
tions as possible, and (ii) to require as little extra computational 
resources as possible. In our model, the emphasis is given to the 
dynamical impact of CRs on hydrodynamics, and not on an ac- 
curate spectral representation of the CRs. The guiding principles 
are energy, pressure, and particle number conservations, as well 
as adiabatic invariants. Non-adiabatic processes wiU be mapped 
onto modifications of these principles. 

The goal of the formalism is to permit explorative studies in- 
vestigating under which circumstances CRs are likely to affect 
galactic or cosmic structure formation processes. The approach 
is phenomenological, which implies that any adopted model pa- 
rameter like CR injection efficiency or diffusivity has to be ver- 
ified by independent means before definite conclusions may be 
drawn, if the conclusions depend strongly on those assumptions. 
Nevertheless, it is plausible that in many cases the overall phys- 
ical picture will not sensitively depend on the details of the CR 
description, but merely on the presence of CRs. An implementa- 
tion of the formalism in a cosinic structure formation code will 
therefore allow an identification of regimes in which CRs are 
likely to have a significant impact on the properties of galaxies 
and galaxy clusters. We want to stress that including "only" an 
approximate description of CR dynamics into galaxy and struc- 
ture formation calculations should be regarded as a significant 
improvement compared to the current situation where the CRs 
are ignored completely. 

In this paper, we give the theoretical basis and justification 
of our formalism. An implementation and first applications of 
our formalism are described by two companion papers. Jubelgas 
et al. (2006) describe the implementation into the GADGET 
code and present first results on galaxy formation including CRs. 
Pfrommer et al. (2006) develop a SPH Mach number measuring 
scheme, which is a necessary prerequisite to study CRs injected 
at accretion shock waves in the cosmic large-scale structure. 

1.2. Approximations and assumptions 

A number of simplifying assumptions are necessary in order to 
reduce the large number of degrees of freedom of a CR popula- 
tion before it can be included in a numerical tractable scheme. 
For transparency, we list the main assumptions here: 

1 . Protons dominate the CR population: We only model the 
dominant CR proton population, assuming that the presence 
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Fig. 1. Thermal and power-law (CR) proton spectra for a variety of tem- 
peratures and spectral indices a. No low momentum cutoff was imposed 
onto the CR populations. Displayed is the energy per logarithmic mo- 
mentum interval over dimensionless proton momenta p = PpKnip c^). 
This form of presentation permits to read-off the dynamically important 
momentum ranges of the populations. All spectra are normalised to ex- 
hibit the same pressure. The thick lines represents a possible location in 
the hot medium of our Galaxy (kT = 10 eV, a = 2.75) where thermal 
equals CR pressure. The break in the CR energy distributions is solely 
due to the change in the kinetic energy formulae in the trans-relativistic 
regime. 



of or-particles and heavier ions would not change the dynam- 
ical picture. In contrast to heavier ions, a-particles carry a 
significant fraction of the total CR energy. Nevertheless, our 
assumption is a reasonable approximation, since the energy 
density of a-particles can be absorbed into the proton spec- 
trum. As a coarse approximation, a GeV energy a-particle 
can be regarded as an ensemble of four individual nucleons 
travehng together due to the relatively weak MeV nuclear 
binding energies compared to the kinetic energy of relativis- 
tic protons. For hadronic interactions, the fact that the four 
nucleons are bound is of minor importance. Since Coulomb 
cooling and ionisation losses are proportional to the square 
of the nucleus' charge, each of the four nucleons of the a- 
particle is experiencing a loss of kinetic energy which is 
identical to the loss that a free CR proton with exactly the 
same specific energy would experience. Furthermore, due 
to their low galactic energy density, but also due to their 
more complex radiative losses, we ignore CR electron pop- 
ulations. Only a description of quasi-equilibrium secondary 
electrons resulting from hadronic CR-gas interactions is pro- 
vided, which has some restricted applicability in galaxy clus- 
ter physics. 

2. A momentum power-law is a typical spectrum: We as- 
sume that the CR momentum spectrum can be well described 
by a power-law with spectral index a. This is not only con- 
sistent with the observation of Galactic CRs (a ^ 2.75) 
but is also predicted by many CR acceleration and diffu- 
sion models, which frequently give momentum power-law 
spectra (see Schlickeiser 2002, for a review). Especially CR 
injection at shock waves in the test-particle limit predicts 
power-law spectra in momentum space. The dynamically rel- 
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evaiit physical quantities of the CR population are its kinetic 
energy density scr, the average energy Tqr = scr/icr. and 
the pressure Pcr. For an assumed power-law spectrum with 
steep spectral index a at high CR momenta, these two quan- 
tities are completely determined by a specification of the nor- 
malisation constant C and the low-momentum cutofl" q of 
the power-law. No high-momentum cutoff of the spectrum 
is considered, since for a sufficiently steep spectrum (a > 2), 
the high-energy range is dynamically unimportant. Instead, 
the dynamics is dominated by particles with momenta clos- 
est to ^ OTp c (the particles at the lower cutoff), or around nip c, 
whichever is larger (and the latter only if a < 3, otherwise 
the particles at q nip c dominate always). 
3. Energy, particle number, and pressure are the relevant 
CR variables: The many degrees of freedom a CR spec- 
trum have to be reduced to simpUfied evolution equations 
for the spectral normalisation constant C, for the cutoff q, 
and eventually for the spectral index a. The mapping of the 
full spectral dynamics on this restricted set of variables is 
not unique. However, physical intuition tells us that the most 
relevant quantities, which should be reproduced in the for- 
malism with highest accuracy, are the CR energy density, 
the number density, and the pressure. CR pressure, which 
is very similar to the CR energy density, will only be used 
if a varying spectral index is chosen. Therefore, whenever 
a physical effect modifies a CR spectrum, the changes on 
C, q, and eventually a will be calculated by taking energy, 
particle and eventually pressure conservation into account. 
Whenever we have to make a choice we will favour accu- 
racy in CR energy over accuracy in CR number density. This 
choice is motivated by our aim to study the dynamical ef- 
fects of a CR population on their host environment. Spectral 
details of the CR populations, which are not the focus of our 
formalism, are therefore treated very coarsely. 

We present two models for the description of the CR popula- 
tion that differ in the degree of their complexity. In the simpler 
model, we assume the CR spectral index a to be the same ev- 
erywhere, while in the more complex model we allow it to vary 
in space and time. The model with the constant spectral index 
suppresses a number of CR phenomena like spectral steepening 
due to transport effects, or spectral flattening due to fresh particle 
injection. However, these effects are probably not essential for a 
first order description of the global dynamical impact of CRs. 

1.3. Captured physics 

Our framework is set up to accommodate a number of essential 
physical processes of CR gases, like particle acceleration, diffu- 
sion, and particle interactions. The different components of our 
framework are of different relevance in the various astrophysical 
environments: 

1. CR energy is an important energy storage within galactic 
metabolisms. In our Galaxy, the CR energy density is com- 
parable to the thermal one. 

2. CR pressure is suspected to play an important role in forming 
galactic winds. 

3. Adiabatic energy losses and gains of CRs are the mechanical 
way to exchange energy with the thermal gas, and therefore 
relevant for CR driven winds or radio bubbles in cluster at- 
mospheres. 

4. Shock acceleration is believed to be an important source of 
CRs. In the galactic context, supernova shock waves should 



dominate the energy injection, in galaxy clusters, the cluster 
merger and accretion shock waves may be more efficient. 

5. The role of CR diffusion is actually not always clear. 
Diffusive escape times in galaxy clusters exceed the Hubble 
time for GeV particles, whereas we know that CRs escape 
from our own galaxy on time-scales of lO' years, which cer- 
tainly involves diffusive CR transport. 

6. In situ CR re-acceleration by plasma waves may energise 
relativistic electrons in several astrophysical environments. 
This process may be important for cluster radio halos but 
also in several other astrophysical situations (Galaxy, Sun,..). 

7. Coulomb and ionisation energy losses of the CRs are heating 
the gas, thereby possibly maintaining a minimum tempera- 
ture and ionisation in dense molecular clouds. 

8. Bremsstrahlung from CR protons is not an important effect 
- to our knowledge -, neither for the CR cooUng, nor for the 
production of high-energy photons. 

9. Hadronic interactions is the only efficient non-adiabatic en- 
ergy loss process of ultra-relativistic protons. Their secon- 
daries (photons, electrons) are important tracers of CR pop- 
ulations. Gamma rays from such interactions are observed 
from our Galaxy. In galaxy clusters, synchrotron emission of 
secondary electrons is proposed as a possible origin of the 
giant radio halos (among other possibiUties, see point 6). 

Thus, these effects should allow realistic studies of the im- 
pact of a variety of physical processes of CRs on galaxies, clus- 
ters of galaxies, and on the large-scale intergalactic medium, in- 
cluding: 

1 . hydrodynamical effects of CRs, 

2. CR injection by diffusive shock acceleration, 

3. in-situ re-acceleration by plasma waves, 

4. non-local feedback from CR injection due to CR diffusion, 

5. CR modified shock structures, 

6. heating of cold gas by CRs, 

7. CR driven galactic winds, 

8. Parker instabiUties of spiral galaxy disks, 

9. rough morphology of gamma ray emission, 

10. rough morphology of radio emission due to secondary elec- 
trons. 

We like to point out that our approach is a phenomenological one 
which does not provide first principle descriptions of the above 
physical effects. It is also not always accurate, especially when 
spectral effects become relevant and deviations from power-law 
spectra are expected, e.g. in the case of non-linear shock ac- 
celeration. Although our formalism captures many properties 
only approximately, it allows to study many astrophysical con- 
sequences in greater detail than before and also permits to study 
the influence of the adopted CR parameters like injection effi- 
ciency, diffusivity, etc. on observable quantities. 

1.4. Structure of paper 

Our basic formalism is outlined in Sect. 2, in which the approx- 
imative description of the CR gas is introduced and its adia- 
batic evolution is described (Sect. 2.1). A model with a con- 
stant spectral index for the power-law description of the CR 
population is given in Sect. 2.2 while the more general ap- 
proach with a spatially and temporally varying spectral index 
is presented in Sect. 2.3. Sect. 3 contains the technical descrip- 
tion of the various non-adiabatic processes: CR injection via 
shocks waves of structure formation (Sect. 3.1.1), and of super- 
nova remnants (Sect. 3.1.2); CR spatial diffusion (Sect. 3.2); in- 
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situ re-acceleration of CRs (Sect. 3.3); CR energy losses due to 
Coulomb interactions (Sect. 3.4), ionisation losses (Sect. 3.5), 
Bremsstrahlung (Sect. 3.6), and hadronic interactions with the 
background gas (Sect. 3.7), including the associated y-radiation 
from the 7r°-decays (Sect. 3.8) and the radio emission of the elec- 
trons and positrons resulting from ;r*-decays (Sect. 3.9). Steady 
state CR spectra, for which injection and cooUng balance each 
other, and freely cooling CR spectra are calculated exactly and 
with our formalism are compared in Sect. 4. Section 5 describes 
how the formalism can be included into a simulation code based 
on smoothed particle hydrodynamics (SPH). Finally, we give our 
conclusions in Sect. 6. 



2. Formalism 

In this section, we develop a description of a CR population 
in a volume element which is comoving with the background 
fluid. This Lagrangian perspective results in a significant simpU- 
fication, since the advective transport processes are fully char- 
acterised by a description of the eff'ects of adiabatic volume 
changes. We will introduce convenient, adiabatically invariant 
variables which are constant in time in the absence of non- 
adiabatic processes. Whennon-adiabatic processes are included, 
they can be expressed in terms of evolutionary equations for the 
adiabatically invariant variables. The chosen formalism is well 
suited for the implementation in Lagrangian numerical simula- 
tion codes, as will be discussed in more detail in Sect. 5 for the 
example of SPH codes. 



2.1. Basic cosmic ray variables 

Since we only consider CR protons, which are at least in our 
Galaxy the dominant CR species, it is convenient to introduce 
the dimensionless momentum p = Pp/(mpc). We assume that 
the difi^erential particle momentum spectrum per volume element 
can be approximated by a single power-law above the minimum 
momentum q: 



dN 
dpdV 



Cp-"eip-q), 



(1) 



where 6ix) denotes the Heaviside step function. Note that we 
use an effective one-dimensional distribution function f(p) = 
A-np^f^^\p). The differential CR spectrum can vary spatially 
and temporally (although for brevity we suppress this in our 
notation) through the spatial dependence of the normalisation 
C = Cix,t) and the cutoff q = q{x,t). In the simpler of our 
two models, we assume the CR spectral index a to be constant 
in space and time, while the more complex model allows for a 
spatial and temporal variation of a as well. 

Adiabatic compression or expansion leaves the phase-space 
density of the CR population invariant, leading to a momentum 
shift according to po ^ p - (p/poy^^ po for a change in density 
from Po to p. Since this is fully reversible, it is useful to introduce 
an invariant cutoff and a normalisation Co which describe the 
CR population via Eqn. (1) if the gas is adiabatically compressed 
or expanded relative to the reference density po- The actual pa- 
rameters of the spectrum are then given by 



q(p) = (p/Po)^ ^0 and C(p) = (p/po) ' Cq. 



(2) 



The variables qo and Cq are hence a convenient choice for a 
Lagrangian description of the ISM or the intra-cluster medium 
(ICM). 



The CR number density is given by 

/-*QO 

ncR= Idpfip) 
Jo 



Cq'-" ^ Coql'" p 
a-l 



a - I Po ' 



(3) 



provided a > I. The kinetic energy density of the CR population 

is 

C rru 



r°° Cm„c 

SCR = \dpf{p)T^{p)= ^ 

Jo a-l 



1„ ( a — 2 3 — a 



\+q^-\ 



(4) 



where Tp{p) = ( Vl + /^^ - 1) nip is the kinetic energy of a 
proton with momentum p. !Bx(a, b) denotes the incomplete Beta- 
function, and a > 2 is assumed. The average CR kinetic energy 
TcR = ecR/ncR is therefore 



■'B 



[ CR 



The CR pressure is 

PCR = ^ jdpfip)pp . 



a — 2 3 — a\ f 



1+q^-l 



S 



a — 2 3 — a] 



(5) 



(6) 



where /3 - v/c - p/ p^ is the dimensionless velocity of 
the CR particle. The CR population can hydrodynamically be 
described by an isotropic pressure component as long as the CRs 
are coupled to the thermal gas by small scale chaotic magnetic 
fields. Note, that for 2 < a < 3 the kinetic energy density and 
pressure of the CR populations are well defined for the limit q 
0, although the total CR number density diverges. 

The local adiabatic exponent of the CR population is defined 

by 



rcR = 



dlogf'cR 



dlogp 



(7) 



where the derivative has to be taken at constant entropy S . Using 
Eqns. (2) and (6), we obtain for the CR adiabatic exponent 



rcR 



P 

PcR [ dC 
a + 2 2 



dPcR dC ^ dPcR dq 



dp dq dp 



I a -2 3 - a\ 



(8) 



Note that the CR adiabatic exponent is time dependent due to its 
dependence on the lower cutoff" of the CR population, q, which 
is shown in Fig. 2. For a mixture of thermal and CR gas, it is 
appropriate to define an eff'ective adiabatic index by 



Tefif = 



dl0g(P,h + /'cR) 



dlogp 



yth -Pth + ycR -Pgr 

^th + ^'CR 



(9) 



2.2. Formalism with a constant spectral index 

The spectrum of the CR population within a given fluid element 
is shaped by a number of different physical processes, such as 
particle injection and escape, continuous and hadronic energy 
losses, or re-acceleration. While all of these processes leave a 
characteristic signature in the CR spectrum, in the framework of 
our simplified model we have to describe their eff'ects in terms of 
the two dynamical variables, C and q (or Co and qo). In order to 
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Fig. 2. Adiabatic index ycR of a CR population as a function of the 
lower momentum cutoff q for a variety of spectral indices a (solid lines). 
Furthermore, the pseudo-adiabatic index Yes. = 1 +^cr/£cr relating the 
total CR pressure and energy density is also displayed (dotted lines). 
Only for a polytropic gas behaviour we have yea. = Tcr- '^^^^ applies 
approximately for CR populations which are either completely non- 
or completely ultra-relativistic. However, in the here important case of 
trans-relativistic CR populations polytropic equations of state can not 
be assumed. 



approximate the key features of the real CR dynamics with our 
description, we have to make a proper choice for how to modify 
our variables by the different processes. The guiding principle 
for this are energy and particle number conservation. 

Consider a non-adiabatic process leading to a change dncR 
in the number density and a change decR in the energy density of 
the particles during an infinitesimal time interval dt. The implied 
change in (C, q) is then given by 



' dC \ / dC/dncR dC/dEcR \ 
dq j \ dq/dncR dq/dscR j 



dncR 
decR 



(10) 



The inverse of the Jacobian can be straight forwardly calculated 
using the definitions in Eqns. (3) and (4) which can easily be 
inverted to yield 



dC 



dq 



, d£cR - Tp(q) dncR 

ecR - Tpiq) ncR 
q decR - TcR dncR 



or - 1 ecR - Tpiq) ncR 



(11) 
(12) 



These relations are reasonable, which can also be verified by the 
following gedanken experiments: if a process increases ecR and 
ncR simultaneously by the same factor 1 + (5, so that decR/scR = 
d«cR/«CR = <5, one gets only a change in the normalisation 
(dC = 6 C), but not in the cutoff (d^ = 0), as it should be. If one 
adds an infinitesimal amount of particles dncR with exactly the 
kinetic energy of the cutoff Tp(q), so that decR = Tp(q) d«cR, the 
normalisation is unchanged (dC = 0), but the cutoff is lowered 
(d^ = -q dncR/[(Q' - 1) "cr] => «cr °^ q^'"), again as expected. 
The adiabatically invariant variables change according to 



\Po 



or+2 
3 



dC = Co 



dscR - Tpiq) dncR 
ecR-Tpiq)ncR 



(13) 



dqo-^ 



dq - 



qo dEcR - TcR d«cR 
a - 1 ecR - Tp(q) hcr 



(14) 



where we used a notation which is mixed in the variant and in- 
variant variables, for convenience. Ways to numerically imple- 
ment the evolution of Co and qo are discussed in Appendix A.l. 

2.3. Formalism with a variable spectral index 

The modehng of certain effects of CR physics such as spectral 
steepening due to transport effects, or spectral flattening due to 
fresh particle injection requires a more elaborate description of 
the CR population. Thus, we generalise the above formalism by 
allowing the spectral index to vary in space and time as well. 
To this end, we now consider the CR population to be described 
by three dynamical variables C, q, and a. In addition to energy 
and particle number conservation, we invoke pressure conserva- 
tion. As a word of caution, we note that the resulting formalism 
is considerably more demanding in terms of computational re- 
sources than the formulation with a constant spectral index. We 
hence suggest that this description is only used if the problem 
under consideration requires this additional information. 

For notational purposes, we define a three-vector with the 
dynamical CR variables as a s (C, q, a)^, where the exponent 
denotes the transpose of the vector. Similarly, we define a vector 
of hydrodynamic CR quantities h = (ncR, ecR, ^'cr)^- Consider 
a non-adiabatic process leading to an infinitesimal change dncR 
in the number density, decR in the energy density, and df cr in 
the pressure of the particles during an infinitesimal time interval 
dt. This implies a change in a = (C, q, a)^ given by 

3 o 

dui = y -^dhi ^ da = Adh. (15) 

7=1 ■' 

Since the entries of the Jacobian A cannot be straightforwardly 
obtained, we propose to compute them by inverting the inverse 
of the Jacobian A"^. Once this is achieved, the changes da in the 
dynamical CR variables can be easily obtained. The inverse of 
the Jacobian is given by 



A-^ = 



ncR 


a- 1 


T 


q 


ecR 


a- 1 


~r 




PCR 


ntp 



ncR 

ncR Tp(q) 
(a -I) 



dncR 
da 
dscR 

da 



p{q)ncR 



CR 



da 



(16) 



where the last colunm can be explicitly expressed as 
dncR 



da 
dscR 
da 

dPcR 

da 



ncR 
a- 1 



[l-l-(a-l)ln^], 



/-*oo 

= - dpTp(p)fip)]np, 
Jo 

= jdppl5{p)f{p)\np. 



(17) 
(18) 

(19) 



We discuss ways to numerically implement the time evolution of 
Co, qo, and a in Appendix A.2. 

A note of caution should be in order. For ultra-relativistic CR 
populations the formalism with varying spectral index becomes 
degenerate, since in this regimes the equation of state is inde- 
pendent of the spectral index (see Fig. 2). However, as can be 
seen in Sect. 4, the dynamics of CR coohng always pushes the 
low-energy cutoff into the trans-relativistic regime, where this 
formalism should work. 
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3. Non-adiabatic processes 

In the following, we discuss various non-adiabatic source and 
sink processes of the CR population. In each subsection, we 
outline the physical motivation of each process and describe 
its effects in terms of a change in energy and number density. 
Afterwards, we generalise to the more complex case of a spa- 
tially and temporally varying CR spectral index. 



3.1. CR shock acceleration 

In this section, we discuss CR acceleration processes at shock 
waves due to gas accretion and galaxy mergers, using the frame- 
work of diffusive shock acceleration. Our description is a mod- 
ification of the approach of Miniati (2001). The shock surface 
separates two regions: the upstream regime defines the region in 
front of the shock which is causally uncoimected for super sonic 
shock waves, whereas the downstream regime Ues in the wake 
of the shock wave. The shock front itself is the region in which 
the mean plasma velocity changes rapidly on small scales, gov- 
erned by plasma physics. In the rest frame of the shock, particles 
are impinging onto the shock surface at a rate per unit area of 
P2V2 - piv\. Here vi and V2 give the plasma velocities (rela- 
tive to the shock's rest frame) in the upstream and downstream 
regimes of the shock, respectively. The corresponding mass den- 
sities are denoted by p\ and p2- 

We assume that after passing though the shock front most 
of the gas thermalises to a MaxweU-Boltzmaim distribution with 
characteristic post-shock temperature T2: 



2nkT2 



ill 



p^ exp 



2 2 

WpC p 

2kT2 



(20) 



where the number density of particles of the thermal distribution 
in the downstream regime, nth = ^2, as well as T2 can be inferred 
by means of the mass, momentum, and energy conservation laws 
at the shock surface for a gas composed of CRs and thermal con- 
stituents. In a companion paper, we describe a formalism for in- 
stantaneously and self-consistently inferring the shock strength 
and all other relevant quantities in the downstream regime of 
the shock within the framework of SPH simulations (Pfrommer 
et al. 2006). Assuming that a fraction of the thermalised particles 
experience stochastic shock acceleration by diff'using back and 
forth over the shock front, the test particle theory of diffusive 
shock acceleration predicts a resulting CR power-law distribu- 
tion in momentum space. Within our model, this CR injection 
mechanism can be treated as an instantaneous process. 

For a particle in the downstream region of the shock to return 
upstream it is necessary to meet two requirements. The particle's 
effective velocity component parallel to the shock normal has 
to be larger than the velocity of the shock wave, and secondly, 
its energy has to be large enough to escape the "trapping" pro- 
cess by Alfven that are generated in the downstream turbulence 
(Malkov & Volk 1995; Malkov & Volk 1998). Thus, only parti- 
cles of the high-energy tail of the distribution are able to return to 
the upstream shock regime in order to become accelerated. The 
complicated detailed physical processes of the specific underly- 
ing acceleration mechanism are conveniently compressed into a 
few parameters (Jones & Kang 1993; Berezhko et al. 1994; Kang 
& Jones 1995), one of which defines the momentum threshold 
for the particles of the thermal distribution to be accelerated. 



^inj — ^injPth — 



2kT2 



m^c 



2 " 



(21) 



Since Coulomb and ionisation losses efficiently modify the low 
energy part of the injected CR spectrum, we propose to fol- 
low the recipe presented at the end of Sect. 3.1.2, i.e. to simply 
increase the low energy spectral break of the actually injected 
spectrum without changing the normalisation of the CR spec- 
trum. 

In the Unear regime of CR acceleration, the thermal distribu- 
tion joins in a smooth manner into the resulting CR power-law 
distribution at (j'inj so that Xj,,] represents the only parameter in 
our simplified diffusive shock acceleration model. 



(22) 



/cR,lin(p) = /th(?mj) I 0(P- qinj)- 

\ ?inj I 



Fixing the normalisation of the injected CR spectrum by this 
continuity condition automatically determines Qnj which de- 
pends on the second adiabatic invariant. The slope of the injected 
CR spectrum is given by 



r+2 



ainj - 



r 



where r = — = — 



Pi 



V2 



(23) 



denotes the shock compression ratio (Bell 1978a,b; Drury 1983). 
In the linear regime, the number density of injected CR particles 

is given by 



Jo 



AncR,lin = I dp/cR,lin(p) = /th(^inj) 



ttinj - 1 ' 



(24) 



This enables us to infer the particle injection efficiency which is 
a measure of the fraction of downstream thermal gas particles 
which experience diffusive shock acceleration. 



'7CR,lin = 



A«CR, 



lin 



inj 



nth 



inj 



1 



■ e ^"i. 



(25) 



The particle injection efficiency is independent of the down- 
stream post-shock temperature T2. These considerations allow 
us to infer the dynamically relevant injected CR energy density 
in the linear regime: 



AecR,im - rjCRMnTcRiamj, qini) nth(T2)- 



(26) 



In our description, the CR energy injection efficiency in the lin- 
ear regime is defined to be the energy density ratio of freshly 
injected CRs to the total dissipated energy density in the down- 
stream regime. 



^lin 



AgCR,lin 
A^diss 



where Aediss = sm - smr''- 



(27) 



The dissipated energy density in the downstream regime, Aejiss, 
is given by the difference of the thermal energy densities in the 
pre- and post-shock regimes, corrected for the contribution of 
the adiabatic part of the energy increase due to the compression 
of the gas over the shock. 

In order to obey energy conservation as well as the equations 
of the linear theory of diffusive shock acceleration, has to 
fulfill a boundary condition which ensures that the dynamical 
pressure exerted by CRs is smaller than the ram pressure pifj of 
the flow, yielding 



^'cR (a - I) c^ T]cRM 



where a = a 



6viV2 



,a-\ 



inj 



a — 2 3 - a 



< 1, (28) 



inj- 
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Fig. 3. CR energy injection efficiency for the diffusive shock accelera- 
tion process. Shown is the CR energy injection efficiency finj (dotted) 
for the three post-shock temperatures kT2/ksW = 0.01, 0.3, and 10. An 
effective CR energy efficiency fcouiomb (solid) is obtained by consider- 
ing Coulomb losses which remove the low-energetic part of the injected 
CR spectrum on a short timescale. In this figure, we adopted the fol- 
lowing values for our model parameters, Xi^j = 3.5, ^max = 0.5, and 

^Coulomb = 0.03. 



In typical applications like cosmological SPH simulations, 
this boundary condition can be substantially simplified. To this 
end, we propose the following modification of the CR energy 
injection efficiency in order to account for the saturation effect 
at high values of the Mach number: 



^inj — 



1 I '1™ 

1 -exp|-- — 

Smax 



(29) 



Numerical studies of shock acceleration suggest a value of 
^max - 0.5 for the limiting case of the CR energy injection ef- 
ficiency (Ryu et al. 2003). One can then infer the injected CR 
energy density in terms of the energy injection efficiency of dif- 
fusive shock acceleration processes. 



AecR,inj = ^injAfidii 



(30) 



We note that nonlinear effects, in the form of back-reactions of 
the accelerated particles upon the shock, change the expectations 
with respect to the linear case. These effects are expected to be 
important well before the limit given by Eqn. (29) (see Eichler 
1979; Drury & Voelk 1981; Axford et al. 1982; Malkov et al. 
2000; Blasi 2002; Blasi et al. 2005; Kang & Jones 2005, and 
many others). 

(DELETED TEXT) 

The average kinetic energy of rcR(Q'inj,^inj) of an injection 
power-law spectrum with CR spectral index amj is given by 
Eqn. (34), but with the lower CR momentum cutoff of Eqn. (21). 
In combination with the slope ainj, the value of Xinj regulates 
the amount of kinetic energy which is transferred to the CRs. 
Theoretical studies of shock acceleration at galactic supernova 
remnants suggest a range of Xinj - 3.3 to 3.6, implying a particle 
injection efficiency of 77cR,iin - 10 to 10"^ (Drury et al. 1989; 
Jones & Kang 1993; Berezhko et al. 1994; Kang & Jones 1995; 
Malkov & Volk 1995). 

Fig. 3 shows the CR energy injection efficiency <finj as a func- 
tion of spectral index ainj . It can be clearly seen that our simpli- 
fied model for the diffusive shock acceleration fails in the limit 



of weak shocks where it over-predicts the injection efficiency. 
Especially in this regime. Coulomb and ionisation losses have to 
be taken into account which remove the low-energetic part of the 
injected CR spectrum efficiently on a short timescale (cf. Sect. 
3.4 & 3.5). This gives rise to an effective CR energy efficiency 
^Coulomb which is obtained by keeping the normalisation of the 
CR spectrum fixed while simultaneously increasing the cutoff: 
^inj — > ^Coulomb-' Thus, in the Unear regime, the effectively in- 
jected CR energy density is given by 

('?C 1 b\'~''^'°' 
.(31) 

The effective CR energy efficiency fcouiomh in the linear regime 
is obtained by analogy to the previous considerations. 



.Coulomb — 



,lin,Coulomb 
^^diss 



(32) 



Following our suggestion for saturation effects of the shock 
acceleration given in Eqn. (29), we can obtain the effec- 
tively injected CR energy density in the non-linear regime 
AecR,inj,Couiomh- Assuming a dominant thermal gas component, 
the spectral index amj can be translated into the Mach number of 
the shock. Mi, depending on the adiabatic index of the thermal 
gas y, 



Ml 



I 2(2-i-ainj) 
1 -I- 2aini ~ 3r 



(33) 



For a pure thermal gas, the spectral index ainj - 2 formally corre- 
sponds to an infinite Mach number. In the case of a varying spec- 
tral index, we re-map the changes AecR.inj, AncR,inj, and the com- 
puted value of the injection spectral index Cinj = (r + 2)/(r - 1) 
onto the dynamical CR variables (C, q, a) which describe the to- 
tal CR population. 

3.1 .1 . CR injection by structure formation sliock waves 

For estimating the CR injection at structure formation shock 
waves in a numerical simulation a dynamical Mach-number 
finder is required, which identifies and characterises shock 
waves on the fly. Such a Mach number finder for SPH simu- 
lations is presented and tested in Pfrommer et al. (2006), and 
applied to the problem of CR injection efficiencies of the differ- 
ent shock waves in the large scale structure. 

3.1 .2. CR injection by supernovae 

Shock waves in supernova remnants are believed to be the 
most important CR injection mechanism in the galactic context. 
However, for typical simulations of galaxy and structure forma- 
tion the spatial scales of supernovae are below the resolution 
length. Therefore we need a sub-grid description for supernova 
CR injection. 

A significant fraction fsN ~ 0.1 - 0.3 of the kinetic energy 
of a supernova may end up in the CR population. Therefore we 
set the energy injection rates into the CR and thermal pools to 
(decR/dOsN = ^SN desN/d? and (deoi/dOsN = (1 - ^sn) desw/d?, 
respectively. Here dssn/dt is the total SN energy release rate 



' Solely for illustration purposes, we sketch how the acceleration ef- 
ficiency is modified for a constant Coulomb cutoft' while we refer to 
Jubelgas et al. (2006) for an algorithm to compute this cutoff on the fly 
in simulations. 
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per volume. The increase in CR number density is given by 
(dncR/d?)sN = (decR/dOsN/T'cR' where 



^ CR 



: OTd C 



Mnj 



or - 2 3 - a 



(34) 



is the average kinetic energy of an injection power-law spec- 
trum with spectral index a\a^ and lower momentum cutoff q^^y 
A plausible value for the injection spectral index is ainj — 2.4 
in a galactic conlext. The low-momentum cutoff can be set to 

^inj ~ .yjkT/imp c^) since the power-law spectrum resulting from 
shock acceleration reaches down to the thermal population with 
temperature kT. 

However, in numerical practice it may be more economical 
to use a somewhat higher value for q-mj, because in many cir- 
cumstances Coulomb and ionisation losses will rapidly remove 
the low energy part of the CR spectrum, so that the energy of 
these CRs is nearly instantaneously reappearing as thermal en- 
ergy shortly after injection. A slight re-calibration of ^sn can 
take this into account, so that a numerical code does not have 
to explicidy follow the appearance of a short-Uved, low energy, 
super-thermal CR population. 

A criterion to find an adequate q^a^ is the requirement that the 
injected CR spectrum above ^inj has a coohng time TcooK^inj, «inj) 
which equals the energy injection timescale Tmj defined as the ra- 
tio of the present CR energy to the energy injection rate (above 
Q'inj).^ The rational behind this approach is that injected CRs with 
P < ^inj would not survive energy losses sufficiently long to 
make a significant change to the existing CR population given 
the current injection rate. Further mathematical details about this 
description can be found in Sec. 4.1.2. 

In the case of a varying spectral index, we also have to cal- 
culate the pressure injection rate Pgr for the injection spectrum 
according to Eq. (6). 

3.2. CR diffusion 

The ubiquitous cosmic magnetic fields prevent charged rela- 
tivistic particles to travel macroscopic distances with their in- 
trinsic velocity close to the speed of light. Instead, the parti- 
cles gyrate around, and travel slowly along magnetic field hnes. 
Occasionally, they get scattered on magnetic irregularities. On 
macroscopic scales, the transport can often be described as a 
diffusion process if the gyro-radius can be regarded to be small. 
The diffusion is highly anisotropic with respect to the direction 
of the local magnetic field, characterised by a parallel k\\ and a 
perpendicular Kj_ diffusion coefficient. Both are usually functions 
of location and particle momentum. 

Microscopically, the scattering of the CR on magnetic irreg- 
ularities of MHD waves slows down the parallel transport, but 
allows the perpendicular transport since it de-places the gyro- 
centre of the CRs. Therefore, both microscopic diffusion coeffi- 
cients depend on the scattering frequency Vscatt> but with inverse 



^ Another criterion could also be given by the requirement that 
the Coulomb-loss timescale Tc(?mj) = |rp(^mj)/(drp(^inj)/dJ)c| (see 
Eqn. (54)) of the particles near the injection cutoff should be on the or- 
der of the simulation time-step. However, this approach has two poten- 
tial shortcomings: First, the dynamics becomes formally time-step de- 
pendent, although the Coulomb cooling dynamics will largely compen- 
sate for this. Second, if small numerical timesteps permit q-,„j q ~ 1, 
the large number of low energy CRs injected may drag an existing CR 
population to lower energies, due to our simplified CR dynamics. 



proportionalities: k\\ oc vl^^^ and k± oc Vscatt- Particles are best 
scattered by MHD waves with a wavelength comparable to the 
CR's gyro-radii, which itself depends on the particle momen- 
tum. The various wavelength bands are populated with different 
strength. Therefore, the scattering frequency is usually a func- 
tion of the particle momentum. The exact functional dependence 
is determined by the plasma turbulence spectrum on scales com- 
parable to the CR gyro-radii. 

In this picture (e.g. see Bieber & Matthaeus 1997), the diffu- 
sion coefficients can be written as 



^Bohm 



e 
s 



I + e 



2 ^Bohm • 



(35) 
(36) 



Here, e «: 1 is the ratio of the scattering frequency Vscatt 
to the gyro-frequency Q. = vlr^, and /^Bohm = ^'''gl^ = 
V pnip c^/(3 Ze B)is the Bohm diffusion coefficient. In most cir- 
cumstances, K± will be many orders of magnitude smaller than 
/C||. Thus, from a microscopic point of view, CR cross field diffu- 
sion seems to be nearly impossible. 

MacroscopicaUy, the cross-field particle transport is much 
faster than the microscopic diffusion coefficient suggests. The 
reason for this is that any small displacement from the initial 
field line, which a CR achieved by a perpendicular microscopic 
diffusion step, can be strongly (often exponentially) amplified if 
the CR travels along its new field line (Rechester & Rosenbluth 
1978; Duffy et al. 1995). This is caused by diverging magnetic 
field lines due to a random walk in a turbulent environment. This 
effect should always be present at some level even if a large-scale 
mean field dominates the general magnetic field orientation. The 
resulting effective diffusion coefficient k± is difficult to estimate 
from first principles (see the discussion in EnBUn 2003), but its 
dependence on the particle momentum is the same as that of 
the parallel diffusion coefficient, due to the dominant role the 
parallel diffusion plays in the effective cross field transport. We 
therefore assume 



K±(p) = 6_lKii(p) 



(37) 



with typically ~ 10 - 10 ^ (Giacalone & Jokipii 1999; 
EnBUn 2003), but see Narayan & Medvedev (2001) for argu- 
ments of a larger 6j_ ~ 10"^ In order to be flexible about the 
underlying MHD turbulence which fixes the momentum depen- 
dency of K\\, we assume 



hip) = hip) = ^oyS/" 7^ = Kop""^' y^-' . 



(38) 



The velocity factor /3 expresses the reduction of diffusion speed 
for non-relativistic particles. For a power-law turbulence spec- 
trum of the form E(k) dk oc "i"* (j/j one obtains dp = 2 - aturb 
and dy = 0, e.g. dp = j for a Kolmogorov-type spectrum. We 
have included the parameter dy in order to allow for low-energy 
deviations from a pure momentum power-law dependence. Such 
deviations can for example be caused by modifications of the 
turbulence spectrum due to MHD-wave-damping by the low- 
energy bulk of the CR population with small gyro-radii. 

It should be noted that turbulent motions in the gas also lead 
to an effective diffusion (e.g. Cho & Lazarian 2004). Due to the 
Lagrangian nature of our CR- fluid description, any turbulent dif- 
fusion due to numerically resolved eddies is automaticaUy in- 
cluded. Unresolved turbulence, however, may have to be taken 
separately into account by adding a momentum independent dif- 
fusion term to Eqn. (38). 
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The equation describing the evolution of the CR spectrum 
f{x, p, t) due to diffusion is 



where the diffusion tensor 



(40) 



is anisotropic with respect to the local main magnetic field direc- 
tion b(x) - B(x)/B{x). Since we are interested in a simplified 
description, we have to translate Eqn. (39) into changes of CR 
number and energy density. Integrating Eqn. (39) over p leads to 
an equation governing the change in CR number density due to 
diffusion: 



I dt 



/difif 



dxi dxj 



C 



2-d„ 



-a+2+d,, 



(1 



+q'y 



I — d-v 



S 



a — I — dp - dy 4 + dp 



(41) 



This can only lead to reasonable results if the condition a > 
1 + dp + dy is fulfilled. For a Kolmogorov-turbulence diffusion 
coefficient {dp = ^, dy - 0), this translates into a > 1.33. 

One could assume that the transported energy is simply 
(decR)diff = TcR (dncR)diff • However, this ansatz would ignore 
that the more energetic particles diffuse faster, implying that the 
effective CR energy diffusion is more rapid than the CR number 
diffusion. In order to model this, we multiply Eqn. (39) by Tp(p) 
and integrate over p. This leads to 



d { C OTp 



dt i^fi dxi dxj\a-2- dp 



-a+2+d„ 



I — dy I a — 1 — dp — dy 4 + dp — a 

— — S_L 

2 i+r 



dy 

2 IT? 



a — 2 — d„ — dy 4 + do — a 



(42) 



This equation can only give reasonable results if the condition 
a > 2 + dp + dy is fulfilled. For a Kolmogorov-turbulence diffu- 
sion coefficient {dp = \, dy = 0), this translates into a > 2.33. 

In the case of a varying spectral index, we need to compute 
{dPcRldt)m which can be obtained by multiplying Eqn. (39) by 
pP{p) and integrating over p. This results in 



X S 



■2-d„ 



■dy 4- 



(43) 



This equation can only lead to reasonable results if again the 



condition a > 2 + dp + dy is fulfilled. 



3.3. CR in-situ re-acceleration 



leading to a re- acceleration of an existing cosmic ray popula- 
tion. Since the CR-number is not changed by this process, we 
can state: 



(39) (dn. 



l dncR \ 



0. 



(44) 



The change in the CR energy can be derived from the Fokker- 
Planck equation for the 3-dimensional momentum distribution 
function p^Hp) - f{p)l{4n p^) of an isotropic CR distribution: 



df\p) 
dt 



j_ d_ 
dp 



p^D 



' dp 



(45) 



where Dp is the (pitch-angle averaged) momentum- space diffu- 
sion coefficient. Dp is a function of p, which we parametrise by 



D„ 



: Do p^'"" y^'"^ ■ 



(46) 



We note that here, as in the section on CR diffusion, a power- 
law spectrum of the relevant plasma waves is assumed, as it is 
expected in the inertial range of MHD turbulence. In case the 
CRs are numerous enough to extract energy from the wave spec- 
trum with a considerable rate (compared to the wave-cascade or 
-decay time of the turbulence), the wave spectrum itself would be 
modified, reducing the acceleration efficiency. A self-consistent 
description would therefore also need to follow ffie wave spec- 
trum and its modification due to the energy extraction by CRs 
(e.g. Miller et al. 1996; Brunetti et al. 2004; Ptuskin et al. 2005; 
Brunetti & Lazarian 2007). However, this is beyond ffie scope of 
ffiis paper. Thus, as soon the cosmic ray population becomes im- 
portant due to re-acceleration, our description has probably left 
its range of validity. 

Taking the appropriate moments of Eqn. (45) leads to evolu- 
tion equations for ffie CR number density, Eqn. (44), and for the 
CR energy density, 

^scr\ ^ {2 + a)CDom„c^ 



dt 



2 W 



a + Op + Gy — 2 2 — a — ap\ 



2 

(l-fl,)/2 



(47) 



which is valid for a > 2 - Up - Uy. This energy has to be taken 
from ffie kinetic-energy dissipation budget of the hydrodynami- 
cal flow. In the case of a varying spectral index, we can obtain ffie 
evolution equation for ffie CR pressure by taking ffie appropriate 
moment of Eqn. (45), 

— — = {2 + a)CDQ—- 

dt /re-acc 3 



2 



a + Qp + Qy 2 — a — Up 



a + Up + Uy — 2 4 — a — Up^ 



+ (^-"-"o (l + q^) 



-a-yil 



(48) 



The diffusive propagation of CRs implies that CR particles scat- 
ter resonantly on plasma waves with wavelength comparable to 
their gyro-radii. Since these waves are propagating, the CRs ex- 
change not only momentum, but also energy with the waves. 



which is valid for a > 2 - - a^. 



The parameterisation of the momentum diffusion coefficient 

is chosen to be similar to the one for the spatial diffusion co- 
efficient because of their related physical background: boffi are 
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due to scattering on the same plasma waves. Quasi-linear cal- 
culations of the Fokker-Planck transport coefficient of charged 
particles interacting with plasma waves (see Schlickeiser 2002) 
demonstrate that both, the spatial- and the momentum-diffusion 
coefficient depend mainly on the pitch-angle diffusion coeffi- 
cient Dftfj = Dfjft{p,fi), where fi = cos9 is the cosine of the 
pitch-angle (9: 



D„ = 



x: 



8 J-, 



2\2 



dju 



(49) 
(50) 



Here, Va is the phase- velocity of the scattering plasma waves 
which are usually assumed to be Alfven waves^ with Va - 
B/ ^J4■np. In the relevant inertia ranges of the MHD-turbulence 
spectra, it is expected ffiat ffie amplitude but not the composi- 
tion of waves changes wiffi wavelength."^ Therefore, the pitch- 
angle diffusion coefficient should be separable in p and yu, e.g. 
D^^(p,lj.) = Di(p) D2(p). This allows us to relate ffie two diffu- 
sion coefficients via 

Dp^^p^VlX2, (51) 
where X2 is a constant of order unity, and is formally given by 



X2 



D2(H) 



(52) 



Therefore, in ffie framework of a quasi-linear approximation, the 
parameters describing in-situ re-acceleration and diffusion are 
related by 



2p — dp , 



■ dy, and Do = V1X2/K0 , 



(53) 



e.g. for a Kolmogorov-like spectrum of Alfven waves Up = dp = 



3 , Oy - dy - 0, and X2 ~ 0(1). 



3.4. Coulomb losses 



The energy loss of a proton with y <si mp/nie by Coulomb losses 
in a plasma is given by Gould (1972): 



dt 



In 



2mf.c^/3p 



'pi 



2 



(54) 



Here, Wpi = -s/AnehiJm^ is ffie plasma frequency, and ric is the 
number density of free electrons. We note that in a neutral gas 
the Coulomb losses can coarsely be estimated wiffi ffie same 
formulae, provided tig is taken to be the total electron number 
density (free plus bound). However, atomic charge shielding ef- 
fects lower ffie Coulomb losses significantly. A more accurate 
description of ionisation losses is given in Sect. 3.5. 

In order to obtain the Coulomb energy losses of the CR 
population, one has to integrate Eqn. (54) over the spectrum 



' Also the stronger damped fast magneto-sonic waves are discussed 
as efficient accelerators due to their higher phase velocity (e.g. Eilek 
1979; Cassano & Brunetti 2005). 

This wont be the case anymore as soon some part of the CR pop- 
ulation extract a significant amount of the wave-energy budget from 
the spectral range they can resonate with. In this case the CR and 
wave spectra become interdependent and a decompositon of the pitch- 
angle diffusion coefficient into a momentum and a pitch-angle depen- 
dent function is impossible (e.g. Miller et al. 1996; Brunetti et al. 2004; 
Ptuskin et al. 2005). 



fip). This integration can certainly be performed numerically. 
For fast and efficient applications, an approximative analytical 
expression might be more practical. We derive such an expres- 
sion by replacing the term fip in the Coulomb logarithm with 
its mean value for the given spectrum, which can be written as 
(fip) = 3 Pculim^ icr)-^ The Coulomb energy losses are ffien 



d£CR \ 

At Ic 



m^c 
2 



In 



2m<,&-(fip) 
-1 a-i 



S 



a - I a\ 
2 '~2) 



(55) 



Since Coulomb losses only affect the lower energy part of 
the spectrum and therefore should leave ffie normalisation un- 
aflFected, we propose to set (dncR/dOc - (<iscR/dt)c/Tp(q). 
This mimics the effect of Coulomb losses on a spectrum quiet 
well, since Coulomb losses remove the lowest energy particles 
most efficiently from the spectrum, moving ffiem into the ther- 
mal pool. Because also their energy is thermalised, we have 

(deth/df)c = -(decR/dOc and (dnp/df)c = -(dncR/d/)c- The 
second equation (which expresses proton number conservation) 
can be neglected for convenience due to the smallness of ffie 
effect. The case of a varying spectral index can be treated in 
the same way as for a constant spectral index because Coulomb 
losses do not significantly modify the spectral index of the ultra- 
relativistic CR population. 

The Coulomb loss time scales r = ecR/(£cR)c of CR popu- 
lations are shown in 4. 



3.5. Ionisation losses 

Ionisation losses of a proton can be calculated wiffi ffie Beffie- 
Bloch equation (Groom & Klein 2000), which we cast into ffie 
form 



dt 



. (56) 



Here, nz is the number density of atomic species with elec- 
tron number Z, and Iz its ionisation potential. b{y) = 



1 + 2yme/mp + (me/mp)^ is a minor correction factor for the 

here relevant regime 7 <sc mp/we, which we ignore in the follow- 
ing. The density correction factor 6z is usually not significant for 
gases, but should be given here for completeness: 



(2y-Dz y\,z<y 
6z^\ 2y-Dz+ azi(yi,z - y)l In 10)*== yo,z <y < yi,z 
I y< yo,z 



(57) 



Here, y = Inp, = I - 2hi(Sa»pi//z) and yo,z, y\,z, «z, kz 
are empirical constants, which characterise the atomic species 
(Sternheimer 1952). The values for Hydrogen and Helium can 
be found in Table 1. It is apparent that the density effect can be 
neglected below 60 GeV, and ffierefore we propose to ignore it 
in applications in which the CR spectrum is trans-relativistic. 

Note ffiat the high energy limit of Eq. (56) is similar, but not 
exactly identical to ffie Coulomb loss formula Eq. (54). 



' This leads to a slight overestimate of the energy losses. A slight un- 
derestimate results by setting (fip) I yj\ + in Eqn. 55. As long 
these two terms lead to similar loss rates our approximately treatment 
is a good description. Otherwise the integration has to be performed 
numerically. 



EnBlin, et al.: Cosmic ray physics in cosmological structure formation 



11 



Element 


Z 


h 




>'i.z 


az 


kz 


H2 


1 


13.6 eV 


1.76 In 10 


3 In 10 


0.34 


5.01 


He 


2 


24.6 eV 


2.0 In 10 


3 In 10 


0.98 


4.11 



1 000.0 F 



Table 1. Atomic parameters as given by Sternheimer (1952). The 
Hydrogen measurements were done with molecular hydrogen, but no 
significant changes for atomic hydrogen are expected with in the accu- 
racy required for our description. 



Similar to the case of the Coulomb losses, we insert char- 
acteristic momenta into the logarithmic factors in order to al- 
low for an analytical integration of the energy losses over the 
CR spectrum. Since the losses are dominated by the trans- 
relativistic regime, we propose to use simply {p^) a (fip) = 
3 Pcr/C^p hcr) and to set j = j In(p^) in order to estimate 6z. 
This yields 



[ dgcR \ 



S 



z 

a - \ a\ 
2 '"2) 



■S 



5z 
2 

a - 



X (58) 



A by-product of the energy losses is the generation of free 
electrons. The production rate of free electrons in a Hydrogen 
gas can be estimated using the empirical observation that per 36 
eV lost by the proton on average a free electron is produced, 
either as a direct or as a secondary knock-on electron (Longair 
1992). Therefore, the ionisation rate is 



1 (dscR\ 
'36eY\ dt jj ■ 



(59) 



The remaining fraction (1 - /i /36eV = 62 %) of the energy loss 
is heating the medium directly. This electron production term 
can be included in any ionisation equilibrium calculation for the 
medium. 



3.6. Bremsstrahlung losses 

Bremsstrahlung energy losses of protons are usually negligible, 
since they are suppressed by a factor ml /nip x 3 10"^ compared 
with electron bremsstrahlung losses, which are in turn usually 
already small. Therefore, we do not include bremsstrahlung en- 
ergy losses of protons in our description. 



3.7. Hadronic losses 

Another important process is the inelastic reaction of CR nu- 
clei with atoms and molecules of interstellar and intergalac- 
tic matter The CR protons interact hadronically with the am- 
bient thermal gas and produce mainly neutral and charged pi- 
ous, provided their momentum exceeds the kinematic threshold 
^thrWpC^ = 0.78 GeV for the reaction. The neutral pions decay 
after a mean lifetime of 9 x 10"'^ s into y-rays while the charged 
pions decay into secondary electrons (and neutrinos): 

TT* ^ yU^" -I- Vft/V/j ^ e* -I- Ve/Ve + + 

With hadronic interactions, only the CR population above the 
kinematic threshold ^thr is visible through its decay products 
in y-rays and synchrotron emission. Because of baryon number 
conservation in strong and electro-weak interactions, we always 
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0.01 
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" a = 2.75 


,.• ■■ yy'a = 2.75 _ 




g = 2.5(r ------ 
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Fig. 4. Cooling times of CR spectra with spectral indices a = 2.25, 
2.50, and 2.75 as a function of the cutoff q (thin solid lines) in a pure 
hydrogen plasma with ne = 1 cm"^. Coulomb and hadronic energy loss 
times are also displayed with dotted and dashed lines, respectively. The 
thick solid lines gives the total cooling time of a mono-energetic CR 
population with p = q. The feature atq - 0.83 is due to the onset of the 
hadronic losses above this momentum. 



end up with pions and two protons in this CR-proton hadronic in- 
teraction (the possibly produced neutron will decay after a mean 
lifetime of 886 s into a second proton). Thus, the CR number 
density is conserved, implying (dncR/df)had - 0. 

The total energy loss of CRs is independent of the detailed 
mechanisms of how the energy has been imparted on pions dur- 
ing hadronic interactions and given by 



d£p , 

— = CnN O-pp^p Tp Oipp - ^thr) . 
°* 'had 



(60) 



Here, o-pp is the total pion cross section which is given by 
Eqn. (64), Kp ^ 1/2 denotes the inelasticity of the reaction in 
the limiting regime (Mannheim & Schlickeiser 1994), and hn = 
«e/(l - jXiif.) denotes the target nucleon density in the ICM, as- 
suming primordial element composition with = 0.24. The 
change in energy density of CRs because of hadronic losses is 
given by 



/had Jo \ df /had 

= cmn d-ppKp SCR (max(^, ^thr)) , 



(61) 



where £cr is given by Eqn. (4) in which the lower spectral break 
q has to be replace by max(q, ^thr)- The case of a varying spectral 
index can be dealt with in the same way as a constant spectral 
index because hadronic losses do not significantly modify the 
CR spectral index. The hadronic loss time scales r = ecR/(£cR)h 
of CR populations are shown in 4. 

3.8. Gamma ray emission 

An analytic formula describing the omnidirectional differential 

y-ray source function resulting from pion-decay of a power-law 
CR population is given in Pfrommer & EnBhn (2004), yielding 
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in our notation: 



SyiEy)dEydV ^ — 



T'C crpp«N 



2Ey 



AEyAV. 



(62) 



The formalism also includes the detailed physical processes at 
the threshold of pion production like the velocity distribution of 
CRs, momentum dependent inelastic CR-proton cross section, 
and kaon decay channels. The shape parameter 6 and the effec- 
tive cross section cTpp depend on the spectral index of the y-ray 
spectrum a according to 

0.14a-^-^ + 0.44 



6 ^ 0.14a-^'' + 0.44 and 
rpp ^ 32-(0.96 + e''''-^''") mbam. 



(63) 
(64) 



There is a detailed discussion in Pfrommer & EnBlin (2004) how 
the y-ray spectral index Uy relates to the spectral index of the par- 
ent CR population a. In Dermer's model, the pion multiplicity 
is independent of energy yielding the relation Uy = a (Dermer 
1986a,b). 

The derivation of the pion-decay induced y-ray source func- 
tion implicitly assumed the kinematic threshold g'thr to be above 
the lower break of the CR spectrum q. This assumption is sat- 
isfied in the case of our Galaxy, where a flattening of the CR 
spectrum occurs below the kinematic energy threshold isthr = 
1.22 GeV (Simpson 1983). If the inequality q < <7thr is violated 
in the simulation for sufficiently long timescales, the resulting 
y-ray emission maps have to be treated with caution. 

Provided the CR population has a power-law spectrum, the 
integrated y-ray source density Ay for pion decay induced y-rays 
can be obtained by integrating the y-ray source function Sy{Ey) 
in Eqn. (62) over an energy interval yielding 



AEy Sy{Ey) 



Ay = Ay{EuE2)= f 

J El 

C m-jfiCCTppm I ntp Y g /f^+l 
ad ntp \2m^j *\ 25 ' 

\ ] 



4( 

3 a6 



a- 1 
IS 



for JG {1,2}. 



(65) 



(66) 



(67) 



3.9. Hadronically induced synchrotron emission 

This section describes the hadronically induced radio syn- 
chrotron emission while employing the steady-state approxima- 
tion for the CR electron spectrum following Dolag & EnBhn 
(2000) and Pfrommer & EnBhn (2004). This is only justified 
if the dynamical and diffusive timescales are long compared to 
the synchrotron timescale. This may well be the case in clusters 
of galaxies, however, probably not in our own Galaxy. Possible 
re-acceleration processes of CR electrons like continuous in- 
situ acceleration via resonant pitch angle scattering by turbulent 
Alfven waves as well as CR electron injection by other processes 
are neglected in this approach. 

Assuming that the parent CR proton population is repre- 
sented by the power-law of Eqn. (1), the CR electron population 
above a GeV is therefore described by a power-law spectrum 



ME,) 
and Ce 



(— y 

VGeV/ 



GeV 

tte - 2 <tt GeV sb + £ph \ GeV 



'pp 



2 4 

mic 



rinC 



(68) 
(69) 



where the efi'ective CR-proton cross section (Tpp is given by 
Eqn. (64), crj is the Thomson cross section, bb - l(%n) is 
the local magnetic field energy density, and Sph = bcmb + estars is 
the energy density of the cosmic microwave background (CMB) 
and starlight photon field, scmb = ^cmb/(^''') expressed 
by an equivalent field strength Bcmb = 3.24 (1 -i- z)^juG. fistars ei- 
ther has to be guessed or calculated from information of the star 
distribution, or ignored. 

The synchrotron emissivity jv at frequency v and per stera- 
dian of such a CR electron population (68), which is located 
in an isotropic distribution of magnetic fields (Eqn. (6.36) in 
Rybicki & Lightman 1979), is obtained after averaging over an 
isotropic distribution of electron pitch angles, yielding 



h 
Be 



A£,yn(«e)Ce 
^^nSB, 

^B, 



SB 



2nml V 



3eGeV^ 



32;r m,c^ I 



r(^) 



(70) 
(71) 

(72) 



where r(fl) denotes the F-function (Abramowitz & Stegun 
1965), ay = (ae - l)/2 = a/2, and B, denotes a (frequency 
dependent) characteristic magnetic field strength which impUes 
a characteristic magnetic energy density sb^. 

4. Testing the accuracy of the formalism 

The rigid spectral form of a single power-law with cutoff, which 
is imposed on the CR populations in our formalism, leads to in- 
accuracies compared to the proper spectral solutions of the evo- 
lution equations of the CR spectra. The question is how large 
are the inaccuracies for different physical quantities, and under 
which circumstances they can or cannot be tolerated. 

Two tests are presented: a steady state situation, where cool- 
ing balances injection and a time dependent, freely cooling CR 
spectrum, without injection. 

4.1. Steady state CR spectrum 

4.1 .1 . Accurate steady state spectrum 

In order to address these questions we analyse a steady state situ- 
ation in which accurate and approximated spectra can be worked 
out analytically. In a homogeneous environment the CR spec- 
trum f(p, f) should follow the evolution equation 



dfip,t) , d fip,t) 
+ T (pip) fiP, t)) = Qip) , 

Ot dp Tiossip) 



where 



Pip) 



dW 
dt 



dW\ 



m' 



(73) 



(74) 



is the momentum loss rate due to Coulomb- and hadronic losses^ 
(see Eqs. (54) & (60)). We assume a fuUy ionized Hydrogen 
plasma with a thermal electron density of ne = 1 cm"^. A time- 
independent CR injection spectrum 



Qip) = e,nj P 



(75) 



* We approximate here the hadronic losses as a continuous energy 
loss. 
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with an injection normalisation Qi^j = 1 Myr~^ (for a here un- 
specified volume element) and ainj = 2.5 is assumed. The en- 
vironment is assumed to be sufficiently extended so that particle 
escape can be ignored, implying T|oss(p) - °o except for the ther- 
malisation of CRs at very low momenta. 
The steady state spectrum is then given by 



m 



QmiP 



l-a-m 



QmjP 



\pip)\ (ainj - 1) (ttinj - l)Ac 



ip' : 
\pl : 



p <^ p* 
p» p* 



(76) 



and its energy distribution can be seen in Fig. 5. The asymptotic 
approximations assume negligible hadronic and Coulomb losses 
in the low and high energy regimes, respectively. Furthermore, 
the weak logarithmic dependence of the Coulomb losses on the 
particle momentum was ignored in the asymptotic equations. 
The cross-over momentum ~ I depends on the ratio of the 
Coulomb to hadronic loss rates coefficients 



In 



3 , 



^had = O-ppKpinpC nN 

and is given by 



Ac Mne^n 

V^had \ We 



eln((2OTec2)/(;iWpl)) 



mp 6-4 cTpp/STp nN 



(77) 
(78) 

a; 1.087. (79) 



Here, we have inserted a fiducial electron density of tie - \ cmT^ 
and assumed a pure hydrogen plasma, p^. ~ 1.1 is accurate to 
10% for electron densities in the range [10-^ 10^] cm ^. As a 
reference case for an analytic approximation to the steady state 
equilibrium CR spectrum we introduce a simple matched asymp- 
totic solution: 



fappmxip} ~ 



QinjP 



(ttinj - l)AcipJ +p-^) ' 

which is also displayed in Fig. 5. 



(80) 



Dividing Eq. (81) by Eq. (83) yields 



Tciq, a) 



TcRiqiaj, ainj) = Tp(q) 1 + 



Thad(^, a) 



(84) 



which is free of any normalisation constant (C, Qnj) and defines 
an impHcit function relating the injection cutoff ^inj and the equi- 
hbrium cutoff q = qiqinj, ainj, a). 

In order to specify q fully a second condition between q and 
^inj is required. A very simplistic, but - as we will see a poste- 
riori - sometimes reasonable choice, is the condition q^nj = q. 
However, our criterion to specify q\nj is the requirement that the 
injected CR spectrum above qi„j has a cooling time Tcooi(?inj, Q'inj) 
which equals the energy injection time-scale Tinj defined as the 
ratio of the present CR energy to the energy injection rate (above 
^inj). Therefore, we have 



Tcooi = Tinj, where 

Tcool ~ T'cool(^inj5 *^inj) 



(85) 



eCR(Cinj,^mj,ffinj) 



Tinj — Tinj (^inj > O^inj 5 C*, ^, Or) 

This leads to the condition 

Tcool(?inj, ttinj) = Tcool(?, a) , 



(d£CR(Cinj, ^inj, C!'inj)/dOc + had 

ecR(C, q, a) 



, and 



(d£CR(Cinj, Q'inj, a;inj)/df)inj 



(86) 



which indeed implies ^inj = ^ in the case ainj = a. In other cases 
one has to solve Eqs. (84) and (86) numerically for qi^j and q. 
The normalisation of the equilibrium CR spectrum can then be 
derived from Eq. (83) and is given by 



C = ginj Tciq, a) 



TcRiq, a) 1 - ainj q 



(87) 



The resulting spectra for ainj = 2.5 and a - 2.25, 2.5, and 2.75 
are displayed in Fig. 5. 



4.1 .2. Approximate steady state spectrum 

In this section, we calculate the equiUbrium spectrum provided 
by our simplified CR formalism. Since our formaUsm is based on 
CR particle and energy conservation, the injected particle num- 
ber and energy should balance the ones due to Coulomb and 
hadronic losses. CR energy balance therefore yields decR/d? = 
(decR/dOinj + (decR/dOcooi = which is equivalent to 



''decR 
i df 



(Q 



> ^inj> fJ^inj) I — 
/inj 



scr(C, q, a) 

Tcooliq, a) 



(81) 



The CR cooling time-scale due to Coulomb and hadronic losses 
is defined by 



had 



ecr(C, q, a) 



(82) 



The injected number of particles must be balanced by the 

number of particles removed from the CR population due to 
Coulomb losses, since hadronic losses conserve particle num- 
ber in the spectrum: dncR/df - (dncR/dOmj + (dncR/dOc = 0. 
Since (dncR/dOc = (decR/dOc/rp(/j) we get 



^ dwcR 



(2inj> ^inj) Q^inj) 



gCR(C, q, a) 
Tp(q)Tc(q,a) 



(83) 



4.1 .3. Comparison of tins spectra 

Although Fig. 5 shows that the approximated and correct CR 
spectra are different in detail, their integral properties are very 
similar. The total energy in the approximated spectrum with 
a - ainj = 2.5 is only 10% above the correct one, and the pres- 
sure difference is even smaller with 5%. Since the low energy 
part of the spectrum is ignored in our formahsm, it is not sur- 
prising that the total CR number density is underestimated by 
60%. Similarly, the Coulomb heating of the thermal background 
gas by the approximated CR spectrum seems to be underesti- 
mated by 80%. However, this is not correct, since in our for- 
mahsm the CRs injected below the cutoff q return their energy 
instantaneously to the thermal gas. If one takes this into account 
one finds that the Coulomb heating rates agree with 1% accu- 
racy. Since our description is energy conserving, this implies that 
the hadronic losses, and thereby the total hadronic induced radi- 
ation in gamma rays, electrons and neutrinos, agree also on a 
1% level. However, the agreement of the differential radiation 
spectra at a given photon, electron or neutrino energy is worse 
(~ 10% . . . 30% overestimate), due to the different normalisation 
of the high energy part of the approximated and the correct CR 
spectra. 

For comparison, we also calculate the accuracy of the 
matched asymptotic spectrum and equilibrium spectra of our for- 
malism for spectral indices different from the injection index. 
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Fig. 5. Spectral energy distribution per logarithmic momentum interval of a steady state CR population within a volume clement of our test system 
in log-linear (left) and log-log (right) representation. CRs arc injected with an spectral index of a-^^j = 2.5. The thick line gives the correct spectrum 
of the full problem, the dashed line is the matched asymptotic solution, and the profiles with the low-momentum cutoff are equilibrium spectra in 
our formalism for the assumed spectral indices a = 2.25 (dotted curve), a = o-inj = 2.5 (solid line), and a = 2.75 (dash-dotted curve). The small 
bump at p = 0.83 in the correct spectrum is due to vanishing hadronic losses below this momentum. The areas under the curves in this plot are 
proportional to the CR energy, which differ less than 20% between approximate and accurate spectra (see Tab. 2). 



Table 2. Relative accuracy of CR energy, pressure, number, and 
Coulomb loss rate in the various spectral approximations displayed in 
Fig. 5, compared to the numerical solution of the problem (also dis- 
played in Fig. 5). 



namely a - 2.25 and a - 2.75. The corresponding relative ac- 
curacies of CR energy, pressure, number, and Coulomb loss rates 
are similar and can be found in Tab. 2. 

This demonstrates that, despite the roughness of the repre- 
sentation of the CR spectrum in our formalism, CR energy den- 
sity and pressure are calculated with an accuracy acceptable for 
first exploration studies of the influence of CRs on galaxy and 
structure formation. The limitations of the approach, especially 
in representing the detailed spectral shape of the CR population, 
have also become clearer. 

4.2. Freely cooling CR spectra 

Setting the right-hand side in Eqn. (73) to zero describes the 
cooling of an initially injected spectrum. Fig. 6 shows the time 
evolution of numerically exact and approximate spectra with an 
initial condition described by (Co,<?o,a) - (1,10^^,2.5) and 
(Co,q(),a) = (1, 10-',2.5), respectively. The numerically exact 
spectra were obtained using the CR spectrum evolution code by 
Jasche et al. (2007, submitted). The approximative CR spec- 
tra were calculated using the implementation of the approxi- 
mative CR evolution by Jubelgas et al. (2006). A fully ionized 
Hydrogen-only plasma with electron density of 1 cm"-* was as- 
sumed. 

For timescales below the cooling time of ~ 50 Myr most 
of the spectral evolution happens in the low-energy cutoff, and 

the approximate treatment captures the evolution of the exact 
solution quite well. Both solution approach a similar asymptotic 



shape, with a low energy cut off close to the trans-relativistic 
regime. However, the cooling rate seems to be different for the 
two descriptions, the approximate solution seems to cool on a 
slightly larger timescale. This inaccuracy is a consequence of 
the approximative description. 

The momentum losses per hadronic interaction are lower 
for trans-relativistic particles than for ultra-relativistic particles. 
In our approximative description the ultra-relativistic particles 
benefit therefore from the reduced energy losses of the trans- 
relativistic CRs. This effect should therefore vanish for an ultra- 
relativistic-only CR spectrum. Indeed, as can be seen in the right 
panel of Fig. 6, as long as all CRs are ultra-relativistic, the ap- 
proximative solution follows very closely the exact one. 

To summarize, within one hadronic cooling time, the approx- 
imative spectral evolution is accurate if the correct spectral index 
is set. However, the long-term evolution of a freely cooling CR 
population becomes inaccurate after a few cooling times. This 
is the more the case the steeper the CR spectrum is. For appli- 
cations where first order effects of the CR dynamics are consid- 
ered, this level of accuracy should be sufficient. For applications 
which require a high level of accuracy over several cooling times 
a proper treatment of the spectral evolution would be needed. 
However, this is beyond the scope of this work. 



5. Smoothed particle hydrodynamics 

In this section, we describe how the dynamical effects of a CR 
population can be included into the smoothed particle hydrody- 
namics (SPH) simulation technique. 
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Fig. 6. Spectral energy distribution per logarithmic momentum interval of a freely cooling CR populations within a volume element of our test 
system. The initial CR populations are described by (Co,?o,«) = (U 10^^,2.5) (left panel) and {Co,qo,a) = (1, 10^,2.5) (right panel). Spectra 
are shown for cooling ages of « (0.01,0.1, 1, 10, 100, 300) Myr. The numerically exact solutions are shown with soUd Unes. The approximative 
solutions are displayed by dotted lines. 



5.1. Lagrangian fluid dynamics 

The Lagrange density of a magneto-hydrodynamical gas-CR 
medium is 

1 , 

£{r, r) = -pr^- sttip. A) - ecR(p, Co, qo) - £b , (88) 

where eth = p^o (p/po)^"' /(t - 1) is the thermal energy density 
of a gas with adiabatic index y and an entropy described by the 
adiabatic invariant A. Any adiabatic invariant X G [Aq, Cq, ^o) is 
simply advected with an adiabatic flow: 



(89) 



dX _ ldX\ 

dt \ d? / non-adiabatic 

where the rhs allows for non-adiabatic changes discussed in 
Sect. 3. The density evolves according to 

dlnp/df = -V-r, (90) 
and the magnetic field according to the MHD-induction law: 



dB 

— = Vx(rxB-?7VxB) 

dt 



5.2. SPH formulation 



(91) 



In smoothed particle hydrodynamics (SPH), the fluid is dis- 
cretised as an ensemble of particles which carry the mass, 
energy, and thermodynamic properties of the fluid elements. 
Macroscopic properties of the medium such as the density at 
position r, of the i-th particle are calculated with adaptive kernel 
estimation in the form 



Pi = Y^mjWi\ri-rj\,hi), 



(92) 



where mj is the mass of the j-th fluid element and W(r, h) is 
the SPH smoothing kernel. The SPH particle positions r, are 
the dynamical variables of the simulation. However, following 



the approach of Springel & Hemquist (2002) which we extend 
in this work to include an additional CR population and to al- 
low for a general equation of state of the gas, the smoothing- 
kemel lengths hi wiU be considered as dynamical variables of a 
Lagrangian function. 

We introduce the CR spectrum of the j'-th SPH particle 



mfiip) = m 



dA^CR m dNcR 



nii 



fiip) 



(93) 



dp dm p(r,) dp dV p{ri) 

with the help of the CR number per momentum and unit gas 
mass fi{p). Our power-law template spectra are then described 
by 



Mp) = Qp-"e(p-qd, 



(94) 



where C, = C(ri)/p(ri) denotes the CR normalisation constant 
of the i-th SPH particle. Similarly, we introduce the CR energy, 
CR density, and thermal energy per unit gas mass with ecR = 
scr/p, "CR = "cr/p, and Sth - Sth/p, respectively. The equations 
defining these quantities and their changes due to adiabatic and 
non-adiabatic processes in terms of C and q can be obtained from 
the corresponding formulae in this paper by replacing C with C. 
For instance, Eqn. (4) yields 



£CR, 



= dpfiip)Tpip) = ^— X 

Jo 



Z 1+qf 



a 

■2 3 - a 



(95) 



and so on. 

The Lagrange formalism provides an elegant way for deriv- 
ing the equations of motions for SPH simulations. The SPH dis- 
cretised Lagrangian is 



(96) 



where e, - eth + ecR is the total energy per mass of the i-th 
SPH particle, and q = ({r,}, [hi], [At]), which should not to be 
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confused with the CR spectral cutoff q, denotes all degrees of 
freedom of the system. These are the components of the SPH 
particle positions r, and the smoothing lengths hi and their ve- 
locities. The /l, s are Lagrange multipliers introduced by Springel 
& Hernquist (2002) in order to incorporate the choice of the 
smoothing length hi into the Lagrangian via the function 



(97) 



where Msph is the required mass within the smoothing kernel. 

Here, we have ignored the description of magnetic fields 
within the SPH-Lagrangian. For the moment, we treat the evo- 
lution of the magnetic field separately from this Lagrangian for- 
malism, adding instead the magnetic forces ad-hoc to the mo- 
mentum equations of the SPH particles. This is along the lines 
of Dolag et al. (1999), and seems to work well in typical cosmo- 
logical settings. However, we note that the dynamical influence 
of magnetic fields can in principle also be included into the SPH 
Lagrange-function as Price & Monaghan (2004a,b) demonstrate. 

If one derives the SPH-equation of motion from a 
Lagrangian, one obtains a dynamical system which obeys energy 
and entropy conservation. Non-adiabatic processes, like shock 
waves, radiative energy losses and energy exchange of the ther- 
mal and relativistic fluids, thermal conduction, and CR diffusion 
have to be added into these equations. The way such process 
should be implemented in case of CR populations should be- 
come clear from this work. 



the SPH particle energy with respect to the local density in terms 
of the total (thermal plus CR) thermodynamical pressure: 



(103) 



One might argue that this derivation should only be correct for 
thermodynamic systems, and therefore not necessarily for CR 
populations which do not exhibit a Boltzmann distribution func- 
tion. However, the concept of entropy, and the concept of adi- 
abatic processes, which do not change phase space density and 
therefore leave entropy constant, is well defined for an arbitrary 
distribution function. Therefore, Eqn. (103) is a generally valid 
result, which can also be confirmed by an expUcit calculation.^ 

Thus, the Lagrange-formalismfor a variable SPH smoothing 
length introduced by Springel & Hernquist (2002) for a poly- 
tropic equation of state can easily be generalised to a general 
equation of state by replacing the thermal gas pressure by the 
total pressure of all fluid components. A calculation along the 
lines of Springel & Hernquist (2002) shows that the SPH-particle 
equations of motion read 



^ = - V 



gi ^V,- Wijihi) + gj ^V,- Wijihj) 



(104) 



with V, = r,. Pi - P,h i + PcR i the thermal plus CR pressure, and 
Wijihi) = Wiln - rjlhi). 



5.3. Equations of motion 

The equations of motions (of the adiabatic, non-magnetic part) 
of the SPH description foUow from the Hamilton principle, 
namely 



=0 

dt dq dq 



(98) 



The equation determining the smoothing length of the j-th 
particle follows from the variation of the action with respect to 
the Lagrange-multiplier /I, . The corresponding part of the Euler- 
Lagrange equations yields 



0, = O, 



(99) 



which for the special form chosen in Eqn. (97) leads to an im- 
plicit equation for hi that has to be solved numerically in prac- 
tice. 

Variation with respect to hi leads to an equation for /I,-: 



(100) 



dei 


d4>i 


dsi dpi 


d^i' 


dhi 


dhi\ 


dpi dhi 


dhi 



Using now Eqn. (97), one gets 
with gi = 



3 m, dsi 
= ——; — g 



1 + 



hi dpi 
3pi dhi 



Anh] dpi 

Furthermore, the thermodynamical pressure is defined as 



P = - 



(dE 

[dV 



(101) 



(102) 



where 5=5, denotes the entropy of a SPH particle volume 

element of size V = V, = p,7ot, and internal energy E = m, e,-. 
This pressure definition can be used to express the derivative of 



6. Conclusion 

We have introduced a simpUfied model for the description of 
cosmic ray physics with the goal of facilitating cosmological hy- 
drodynamical simulations that self-consistently account for the 
dynamical eflfects of CRs during structure formation. We have 
shown how various adiabatic and non-adiabatic processes can be 
described in terms of a simple model for the local CR spectrum, 
consisting of a power-law with varying normalisation and low- 
energy cutoff. In our basic model, the CR spectral index is held 
fixed and chosen in advance to resemble a typical spectral index 
for the system under consideration. We have also described an 
extended model where the spectral index is allowed to vary as 
well, which leads to a numerically more involved scheme. 

We have demonstrated that dynamical quantities like CR en- 
ergy, pressure, and energy loss rates are all reasonable well rep- 
resented by our approximative treatment of the CR spectra in 
case of a steady state situation in which CR injection and cooling 
balance each other. The accuracy is ~ 10% even if the spectral 
indices of CR injection and population do not match. Given the 
large uncertainties in our knowledge of the parameters determin- 
ing the CR physics like CR injection efficiency, the level of accu- 
racy of our numerical treatment seems to be sufficient. Thus, the 
formalism is suitable for explorations of the possible dynamical 
impact of CRs on galaxy and large-scale structure formation. 

We have explained how our description of the CR physics 
can be self-consistently included into the SPH simulation 
methodology. The formulation we propose manifestly conserves 
energy and particle number. In particular, CR entropy is exactly 



' The internal energy per SPH particle of an ideal (thermal and/or rel- 
ativistic) gas can be written as w, Sj = 2o J dpfajip) Ta{p), where 
a is the index over the particle species (electrons, protons, etc.) with 
momentum-space distribution functions f„j(p), and Ta(p) the relativis- 
tic correct Icinetic energy of the particles (Eqn. (5)). A straightforward 
calculation of dsi/dpi, which uses the first equality in Eqn. (6), leads 
then to Eqn. (103). 
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conserved in adiabatic processes, and the dynamical forces from 
CR pressure gradients are derived from a variation principle. 

This paper has outUned a basic framework for future work 
on the impact of CR populations on galaxy and large-scale struc- 
ture formation. It is accompanied by two papers describing the 
implementation and testing of (i) the CR formalism as described 
here in the GADGET simulation code (Jubelgas et al. 2006), and 
(ii) a shock capturing method for SPH that allows on-the-fly esti- 
mation of the Mach number of structure formation shock waves, 
which is essential to follow CR injection accurately (Pfrommer 
et al. 2006). Further science applications are in preparation. 

Appendix A: Numerically updating the CR spectrum 

A.1. Constant spectral index 

Updating the adiabatic invariant variables Co and qo is most con- 
veniently done via Eqn. (13) and (14). However, if the relative 
changes during a numerical time-step are large, e.g. due to rapid 
CR production at a location without a substantial initial CR pop- 
ulation, these equations would have to be integrated on a refined 
time-grid, or solved with an implicit integration scheme. Both 
methods would be very time-consuming. Therefore, we propose 
another updating scheme: from the initial variables C'o(?o) and 
qo(to) at time to, the corresponding instantaneously particle num- 
ber ncR, energy density ecR, and average particle energy Tcr 
are calculated according to Eqns. (2) to (5). Then hcR, scr, and 
Tqr are updated using the non-adiabatic CR energy and number 
losses or gains during that time-step. And finally, these updated 
values have to be translated back into updated values of CoC^i) 
and ^o(fi). This is easiest by first inverting Eqn. (5) in order to 
calculate q, and then to use Eqns. (2) and (3) to get the updated 
Co(/i) and qo(ti). The inversion of Eqn. (5) has to be done nu- 
merically for Tqr ~ ffJp c^, e.g. using pre-calculated numerical 
tables. However, for the asymptotic regimes we propose the fol- 
lowing approximate inversion formulae: 



^(T+1), T = TcR/(mpc2)» 1 
with q. = {-\ ,andS = S^,^ 



(A.l) 



(A.2) 



A.2. Variable spectral index 

If the relative changes of the dynamic CR variables C, q, and a 
during a numerical time-step are large we propose the following 
numerically efficient updating scheme (instead of updating via 
Eqns. (15) and (16)): from the initial variables Co(fo), qoito), and 
a(to) at time fo, the corresponding instantaneous particle num- 
ber ncR, energy density scr, and pressure Pcr are calculated 
according to Eqns. (2) to (6). Then hcR, ecR> and Pcr are up- 
dated according to the losses or gains of non-adiabatic CR en- 
ergy, pressure, and number density during that time-step. And 
finally, these updated values have to be translated back into up- 
dated values of Co(fi), qoOi), and aiti). This can be done by 
numerically solving the following equation for the root q, 



Pcr q'-"^"^ 



«CR aiq) - 1 
a{q)-l 



-S 1 
IT? 

SPCR 



a(q) -2 3 - a(q) 



)• 



ecR - Tpiq)ncR 



and (A.3) 



(A.4) 



These equations are obtained by combining Eqns. (3), (4), and 
(6). The new CR spectral index a and C are obtained by 
Eqns. (A.4) and (3). 
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